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SUMMARY 

A comparative study was conducted for computational fluid dynamic solutions to flow 
over a backward-facing step. This flow is a benchmark problem, with a simple geometry, but 
involves complicated flow physics such as free shear layers, reattaching flow, recirculation, and 
high turbulence intensities. Three Reynolds-averaged Navier-Stokes flow solvers with k-e 
turbulence models were used, each using a different solution algorithm: finite difference, finite 
element, and hybrid finite element - finite difference. Comparisons were made with existing 
experimental data. Results showed that velocity profiles and reattachment lengths were predicted 
reasonably well by all three methods, while the skin friction coefficients were more difficult to 
predict accurately. It was noted that in general, selecting an appropriate solver for each problem to 
be considered is important. 


INTRODUCTION 

In the application of computational fluid dynamic (CFD) methods, some considerations 
are the relative accuracy, efficiency, and applicability of different flow solvers. Selecting an 
appropriate solver for a given problem based on those criteria is an important factor in the 
successful use of CFD. Conversely, selecting an inappropriate solver can result in an inaccurate 
solution, a waste of computer resources, excessive user effort, a solution with insufficient detail, a 
solution with overwhelming detail, or a failure to obtain any solution at all. 

Validation is a first step in applying CFD to any problem. This is accomplished by 
running the solver on appropriate benchmark cases and comparing the solutions with experimental 
results, or better yet, by identifying such validation cases that have already been performed. The 
benchmark cases, which can vary in complexity and size, should encompass the main fluid 
dynamic features contained in the problem to be investigated. 

For air breathing propulsion systems, separated-reattaching flows are of interest, especially 
as they apply to areas such as diffusers and nozzle boattails, because the presence and 
characteristics of separated-reattaching flow on these components can drastically alter their 
performance. In these flows, the separation is driven by boundary layer development and pressure 
gradients in the flow over a complex geometry, so predicting the separation point is a difficult 
t3.sk 

A classic benchmark case used for validating CFD codes for turbulent separated- 
reattaching flows is the flow over a backward-facing step (BFS). The geometry of the 
configuration is simple, but the flow is still complex and challenging to simulate, involving flow 
physics such as free shear layers, reattaching flow, and recirculation. Proper modeling of 
turbulence is needed for accurate simulation of the flow, and many studies applying different 
turbulence models to BFS flows have been conducted, such as in references 1-3. In BFS flow. 
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the separation point is fixed at the comer of the step; the separation is not driven by boundary layer 
development and adverse pressure gradients. Thus the problem is simplified considerably, 
facilitating study of the reattachment itself. 

In order to assess the application of different CFD solvers to separated-reattaching flows, a 
comparative study was conducted for solutions to flow over BFS. Three Reynolds-averaged 
Navier-Stokes flow solvers incorporating two-equation k-£ turbulence models were used, each 
one using a different solution algorithm: a finite difference method (FDM) solver originally 
oriented to aeropropulsion flows, a computer aided engineering (CAE) oriented finite element 
method (FEM) code, and a CAE oriented hybrid finite element - finite difference method (HM) 
solver. FDM used a low Reynolds number k-s turbulence (LR k-£) model, whereas FEM and 
HM used versions of high Reynolds number k-£ turbulence (HR k-£) models. Solutions were 
compared with each other, and with existing experimental data. This is not a comparison between 
numerical algorithms of finite differences and finite elements, due to a number of factors, and 
particularly because of differences in the turbulence models used. 


EXPERIMENTAL STUDY 

The experimental study of flow over BFS conducted by Driver and Seegmiller (4) was 
selected for use in the present work, because of the extensive quantitative measurements made, the 
relatively high flow speeds, and the low streamwise pressure gradient of the flow. The freestream 
velocity was 145 ft/s (Mach number = 0. 128) in standard atmosphere, and the step height was 0.5 
inches, giving a step height Reynolds number of 33420. The Reynolds number of the boundary 
layer momentum thickness 4 step heights upstream of the step was 5000, producing a fully 
turbulent boundary layer. The ratio of step height to tunnel exit height was 1 :9 for the parallel wall 
case, so the streamwise pressure gradient due to the step was relatively small. Mean and turbulent 
flow velocity components were obtained using two component laser doppler velocimetry (LDV), 
wall static pressures downstream of the step were measured by static taps, and wall friction 
coefficients were obtained from oil-flow laser interferometry. Turbulence intensities, correlations, 
and dissipation were computed from LDV data. 


COMPUTATIONAL METHODS 

Three Reynolds-averaged Navier-Stokes flow solvers with k-£ turbulence models were 
used to solve BFS flow, each using a different solution algorithm: FDM with LR k-£, FEM with 
HR k-£, and HM with HR k-£. HR k-£ assumes law-of-the-wall profiles in the boundary layer, 
whereas LR k-£ requires a large number (20+) of grid points near the wall to numerically resolve 
the boundary layer. Therefore, although LR k-£ is much more computationally intensive than HR 
k-£, LR k-£ is usually expected to be more accurate for separated flows, because it does not 
assume the law-of-the-wall profiles that are based on attached flat plate boundary layers. The 
computational domain, schematically depicted in figure 1, begins four step heights upstream of the 
step where an incoming boundary layer velocity profile is specified, and ends far downstream of 
the step. HM and FEM use the same grid; FDM uses a much finer grid in order to accommodate 
the LR k-£ model. 


Finite Difference Method (FDM) 

The finite difference solver used in this study is PARC, an internal flow Navier-Stokes 
code used extensively by government and industry to analyze propulsion flows (5,6). Two- 
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dimensional / axisymmetric (2D), and three-dimensional (3D) versions are available. PARC was 
derived from the ARC Navier-Stokes code (7) which has been used for external flow calculations. 

The governing equations of motion are the time dependent Reynolds averaged Navier- 
Stokes equations using a perfect gas relationship and Fourier’s heat conduction law. These 
equations are discretized in conservation law form with respect to general curvilinear coordinates 
and solved with the Beam and Warming approximate factorization algorithm (8). Although the 
time dependent formulation of the equations is used, the code is intended for steady state flow 
simulations. The LR k-£ model of Chien (9), including modifications for compressibility by 
Nichols (10), was used in these calculations using the 2D version of the code. 

A 11,341 point grid was used, with 111 points in the horizontal direction and 131 points in 
the vertical direction downstream of the step. In the boundary layers, the grid was packed to the 
walls such that downstream of the reattachment point, on the average, the first grid point off the 
wall was approximately at y+ ~ 1.0. 

Finite Element Method (FEM) 

The finite element solver used is FIDAP Version 5.04 by Fluid Dynamics International 
(11). It is intended to handle a wide variety of applications, including: internal and external 
aerodynamics, heat transfer, crystal growth, extrusion, injection molding, chemical mixing, 
chemical vapor deposition, etc. The code includes its own pre and post-processors to form a 
complete, integrated analysis package. Files can also be imported from, or exported to, other finite 

element pre/post-processors such as PATRAN, ANSYS, and SUPERTAB. 

FEM uses the incompressible, Reynolds-averaged Navier-Stokes equations in either 
dimensional or non-dimensional form. Non-inertial frames of reference are supported, as are free 
surfaces. The Galerkin method of weighted residuals is used for finite-element formulation of the 
problem. A streamline upwinding capability is included to handle the numerical instabilities posed 
by Galerkin's method. Also, the standard pressure discretization can be replaced by a penalty 
function approximation. Several methods are available for solving the discrete system of 
equations including successive substitution, Newton-Raphson, and modified Newton-Raphson 

methods. . , . . . , ,, TTO 

The HR k-e model of Launder and Spalding (12), generally considered the standard HR 

k-e model, is used. However the treatment at the wall differs from that used in the HM (to be 
discussed later). In order to solve the full set of elliptic equations down through the viscous 
sublayer to the wall, a one-element thick layer of special wall elements are used adjacent to the 
physical boundary. These wall elements specify shape functions based on the law-of-the-wall 
profiles to characterize the variations in the mean flow variables. Then, van Driest s mixing length 
approach is used to model the variations of turbulent quantities. Alternatively, a user-defined 
algebraic turbulence model (based on mixing-length theory) could be created and used in place of 
HR k-e; however, this option was not exercised in the present study. 

The computational grid was generated using the mesh generator integrated with the FEM 
system. With the exception of the additional wall elements, the grid is identical to the one used 
with HM, consisting of 2751 nodes and 2796 elements, with 70 elements in the horizontal 
direction, and 40 elements in the vertical direction downstream of the step. 

Hybrid Method (HM) 


FLOTRAN version 2.0, by Compuflo, Inc. (13) uses a hybrid finite-element/ finite- 
difference method. The solution algorithm is based on finite element methods. However, certain 
modifications make the computational efficiency and storage requirements more competitive with 
finite difference / finite volume codes, while still retaining the geometric flexibility of finite 
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element codes. This program consists of the flow solver only; grid generation and post 
processing are handled by finite element pre/postprocessors, such as PATRAN, ANSYS and I- 
DEAS. 

The governing equations are the steady, Reynolds averaged Navier-Stokes equations, used 
in the incompressible form for the present problem. Galerkin’s method of weighted residuals is 
used to discretize the diffusion and source terms. A monotone streamline upwind method is used 
to discretize the advection terms, in order to avoid oscillations in the advection terms, called 
dispersion errors, caused by the Galerkin’s method. The intent of streamwise upwinding is to 
minimize the ‘crosswind’ numerical diffusion that may occur when the fluid flows diagonally 
across the cell. Pressure is derived from a velocity-pressure relation. 

As in FEM, the Launder and Spalding HR k-e turbulence model is used (11). Analytical 
law of the wall and log law of the wall velocity profiles are used, but the exact treatment near the 
wall differs from that used in FEM. 

The grid used by HM in the present study was generated using PATRAN. In order to 
check the grid sensitivity, two grids were created: a coarse grid with 2640 elements, and a fine grid 
with about 6000 elements. However, both produced nearly identical solutions, therefore results 
shown herein are with the coarse grid, which is nearly identical to the grid used by FEM. 


RESULTS 
Reattachment Length 

A sensitive parameter often used to quantify the accuracy of solutions to BFS flow is the 
reattachment length, defined as the distance downstream of the step where the flow separating at 
the step comer reattaches with the bottom wall. Reattachment occurs where the velocity gradient 
off the wall is zero, or in other words, where the wall shear stress is zero. However, there appears 
to be an inconsistency in the experimental results (4). Driver and Seegmiller report a reattachment 
length of x/H = 6.25, as measured by laser interferometry. On the other hand, velocity profiles 
measured by LDV indicate a reattachment point just downstream of x/H = 5, and the profile at 
x/H=6 appears clearly attached. Driver & Seegmiller speculate the discrepancy may be due to a 
long, thin separated region along the wall with very small length scales; but no such flow was 
indicated even in FDM results using LR k-e (which is expected to better resolve near-wall flows 
than HR k-e). Based on the reattachment length determined from oil-flow interferometry, FDM 
gives a good prediction, whereas FEM and HM underpredict by about 12% and 17%, respectively 
(table 1). However, if one were to consider the reattachment length inferred by the LDV velocity 
profiles, FEM and HM give the better predictions. 

Wall Shear Stresses 

Present results concur with studies in the literature (1-3) that wall shear stresses 
downstream of the step are difficult to predict accurately, yet in general, accurate prediction of skin 
friction is important for calculating the drag of practical configurations. Figure 2 depicts the 
present results. FDM created significant under and over shoots of skin friction coefficients, 
behind and ahead of the reattachment point, respectively. FEM produced a reasonable looking 
curve, but it was slightly shifted to the left due to a short reattachment length. HM underpredicted 
the skin friction, and the curve was also shifted to the left due to a short reattachment length. 
Neither FEM nor HM appeared to reproduce the secondary recirculation zone near the base of the 
step, indicated experimentally by positive C f values in that region. 
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Mean Velocity Profiles 


Axial velocity profiles were predicted reasonably well by all three methods, at the stations 
x/H = 3, 6 and 12 downstream of the step (fig. 3). At x/H = 3, all three solutions look 
qualitatively similar, with the ‘comer’ in the velocity profile at the correct y height. At x/H = 6 and 
x/H = 12, the velocity profile comers are still at the correct heights, but all three methods slightly 
underpredict the velocity near the wall, with FDM having the greatest discrepancy, and HM the 
least. 


DISCUSSION 

Turbulence is a major phenomena in BFS flow, and as expected, turbulence models have 
been observed to significantly impact the CFD solutions. Because the turbulence models used in 
FDM, FEM and HM are not exactly the same, differences in the solutions are likely due at least in 
part to different turbulence models. Using LR k-£, Avva, et al. (1) observed strong sensitivity of 
the skin friction results to grid refinement in the inner layer, by varying the number of grid points 
in the inner layer between 5 and 30. In the present study, for FDM with LR k-£, approximately 
18 grid points were in the inner layer, so the skin friction results may change with changes in the 
boundary layer grid. Conventional wisdom holds that HR k-£ is not accurate for separated flow, 
because it assumes an attached velocity profile at the wall. Failure of solvers equipped with HR k- 
£ (FEM and HM) to reproduce the secondary recirculation zone at the base of the step seems to 
support this assertion. However, when the overall situation is considered, the present study 
concurs with Awa, et al. (1) and Steffen (3) that high Re k-£ may actually be preferable to LR k-£ 
for some cases such as BFS flow, because HR k-£ uses less computer resources, is less grid- 
resolution sensitive, and produces results of similar quality overall, when compared to LR k-£. 

Computational efficiency is a consideration, since these codes may be applied to large, 
complex problems; to this end, simple benchmark cases should not consume large quantities of 
computer resources. To obtain the present solutions, FDM used 50,000 Cray Y-MP CPU 
seconds, FEM consumed 615 Y-MP CPU seconds, and HM used 2000 seconds on a VAX9000. 
FEM and HM appear to be substantially more CPU efficient than FDM, but some other factors 
need to be considered. This particular FDM code is known to be very slow to converge for low 
speed flows (freestream Mach number < 0.2). Convergence is much faster for higher Mach 
numbers that are more typical of aerospace applications. Furthermore, as discussed above, LR k- 
£, US ed in FDM, is inherently more computer intensive than HR k-£, used in FEM and HM. 
Therefore, if the problem under consideration has higher freestream velocities as is more typical of 
aerospace applications, or if a HR k-£ model is successfully integrated with FDM, the relative 

computational efficiency of FDM should improve substantially. 

The present study is only one test case, with a single fixed flow regime and a certain set of 
flow features; to draw generalized conclusions regarding the merit or performance of the three 
codes based on this one case only would be misleading. The original objective of this work- 
which was not fully realized due to time and code capability constraints-was to investigate a 
number of geometrically simple benchmark cases encompassing most of the flow regimes and 
flow physics encountered in air breathing propulsion inlet and nozzle flows, using a number of 
different flow solvers. Compressible and supersonic flows are aspects of those flows, therefore a 
simple low angle ‘bent wall’ was selected as a benchmark case to evaluate solver performance at 
high subsonic and supersonic Mach numbers. However, the FEM version available at the time 
did not handle compressible or supersonic flow, and attempts to run this case on HM were 
unsuccessful. Therefore, it becomes clear that an appropriate flow solver must be chosen for each 
problem to be investigated, and that there is no one code that is ideal for solving all types of flow. 
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Table 1 . Reattachment lengths 


Method 

reattachment length (x/H) 

Driver-SeegmiUer Data 

6.25 

FDM 

6.26 

FEM 

5.47 

HM 

5.22 
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SECONDARY 

RECIRCULATION 


Figure 1. Schematic diagram of backward-facing step 
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Figure 2. Wall skin friction coefficient distributions 
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Figure 3. Mean velocity profiles in the separated and reattached regions 
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